{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "27f04776-7c20-4486-ab05-0ce2d2d8a6fb",
   "metadata": {},
   "source": [
    "## 逻辑回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b51d6857-cfd2-476a-96f3-f101cf4aeccb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# options(repos = c(CRAN = \"https://mirrors.aliyun.com/CRAN/\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d0f756eb-99db-4089-8593-30f67745a3f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# install.packages(\"AER\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e713efba-9e3b-4f62-80e1-76e7c2fceb0f",
   "metadata": {},
   "source": [
    "jupyter notebook无法加载数据，直接在R软件里面加载数据，并且保存成csv文件后读入数据\n",
    "```r\n",
    "data(Affairs,package=\"AER\")\n",
    "write.csv(Affairs,\"Affairs.csv\",row.names = FALSE)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "49586704-56b2-4f4b-bb3c-43eceea1b061",
   "metadata": {},
   "outputs": [],
   "source": [
    "Affairs <- read.table(\"data/Affairs.csv\",header=TRUE,sep=',')#sep一定要写，不然可能出错"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "537f1e63-043a-46a8-97e1-3103e913fb29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  affairs gender age yearsmarried children religiousness education occupation\n",
      "1       0   male  37        10.00       no             3        18          7\n",
      "2       0 female  27         4.00       no             4        14          6\n",
      "3       0 female  32        15.00      yes             1        12          1\n",
      "4       0   male  57        15.00      yes             5        18          6\n",
      "5       0   male  22         0.75       no             2        17          6\n",
      "6       0 female  32         1.50       no             2        17          5\n",
      "  rating\n",
      "1      4\n",
      "2      4\n",
      "3      4\n",
      "4      5\n",
      "5      3\n",
      "6      5\n"
     ]
    }
   ],
   "source": [
    "print(head(Affairs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fc7b40cf-ca27-4036-b876-a1608edc910f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "    affairs          gender               age         yearsmarried   \n",
       " Min.   : 0.000   Length:601         Min.   :17.50   Min.   : 0.125  \n",
       " 1st Qu.: 0.000   Class :character   1st Qu.:27.00   1st Qu.: 4.000  \n",
       " Median : 0.000   Mode  :character   Median :32.00   Median : 7.000  \n",
       " Mean   : 1.456                      Mean   :32.49   Mean   : 8.178  \n",
       " 3rd Qu.: 0.000                      3rd Qu.:37.00   3rd Qu.:15.000  \n",
       " Max.   :12.000                      Max.   :57.00   Max.   :15.000  \n",
       "   children         religiousness     education       occupation   \n",
       " Length:601         Min.   :1.000   Min.   : 9.00   Min.   :1.000  \n",
       " Class :character   1st Qu.:2.000   1st Qu.:14.00   1st Qu.:3.000  \n",
       " Mode  :character   Median :3.000   Median :16.00   Median :5.000  \n",
       "                    Mean   :3.116   Mean   :16.17   Mean   :4.195  \n",
       "                    3rd Qu.:4.000   3rd Qu.:18.00   3rd Qu.:6.000  \n",
       "                    Max.   :5.000   Max.   :20.00   Max.   :7.000  \n",
       "     rating     \n",
       " Min.   :1.000  \n",
       " 1st Qu.:3.000  \n",
       " Median :4.000  \n",
       " Mean   :3.932  \n",
       " 3rd Qu.:5.000  \n",
       " Max.   :5.000  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "summary(Affairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "57bd12a0-baa6-4874-a91e-ae18982eca99",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  0   1   2   3   7  12 \n",
       "451  34  17  19  42  38 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table(Affairs$affairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fcfe600e-72bb-45bf-93f8-bfef29dc0d89",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "female   male \n",
       "   315    286 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table(Affairs$gender)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1e3890ba-f00d-4169-8c37-9ebc1fd919dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  gender  affairs      age\n",
      "1 female 1.419048 30.80159\n",
      "2   male 1.496503 34.34441\n"
     ]
    }
   ],
   "source": [
    "print(aggregate(Affairs[c(\"affairs\",\"age\")],by=list(gender=Affairs$gender),FUN=mean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "03ee772f-35e0-4b55-bf5a-b7c7caee471f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " No Yes \n",
       "451 150 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Affairs$ynaffair[Affairs$affairs>0] <- 1\n",
    "Affairs$ynaffair[Affairs$affairs==0] <- 0\n",
    "Affairs$ynaffair <- factor(Affairs$ynaffair,levels=c(0,1),labels=c(\"No\",'Yes'))\n",
    "table(Affairs$ynaffair)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "401c7130-43c1-43f3-8a5b-c444b80a0dd1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = ynaffair ~ gender + age + yearsmarried + children + \n",
       "    religiousness + education + occupation + rating, family = binomial(), \n",
       "    data = Affairs)\n",
       "\n",
       "Coefficients:\n",
       "              Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)    1.37726    0.88776   1.551 0.120807    \n",
       "gendermale     0.28029    0.23909   1.172 0.241083    \n",
       "age           -0.04426    0.01825  -2.425 0.015301 *  \n",
       "yearsmarried   0.09477    0.03221   2.942 0.003262 ** \n",
       "childrenyes    0.39767    0.29151   1.364 0.172508    \n",
       "religiousness -0.32472    0.08975  -3.618 0.000297 ***\n",
       "education      0.02105    0.05051   0.417 0.676851    \n",
       "occupation     0.03092    0.07178   0.431 0.666630    \n",
       "rating        -0.46845    0.09091  -5.153 2.56e-07 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 675.38  on 600  degrees of freedom\n",
       "Residual deviance: 609.51  on 592  degrees of freedom\n",
       "AIC: 627.51\n",
       "\n",
       "Number of Fisher Scoring iterations: 4\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children + religiousness + education + occupation + rating,\n",
    "                data=Affairs,family=binomial())\n",
    "summary(fit.full)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c4ebdcb-e03a-412d-bcc7-bb22bf83483e",
   "metadata": {},
   "source": [
    "残差偏差(Residual deviance)\n",
    "\n",
    "Residual Deviance 表示当前模型与饱和模型（Saturated Model）之间的对数似然差异，公式为：$\\text{Residual Deviance} = -2 \\left( \\ell(\\hat{\\boldsymbol{\\theta}}) - \\ell(\\hat{\\boldsymbol{\\theta}}_{\\text{sat}}) \\right)$\n",
    "\n",
    "其中：\n",
    "- $\\ell(\\hat{\\boldsymbol{\\theta}})$ 是当前模型的对数似然函数值（参数为模型估计值），\n",
    "- $\\ell(\\hat{\\boldsymbol{\\theta}}_{\\text{sat}})$ 是饱和模型的对数似然函数值（参数使模型完美拟合数据）。\n",
    "\n",
    "本质：Deviance 是广义线性模型中 “残差平方和” 的推广，值越小表示模型拟合越好。\n",
    "\n",
    "假设数据 $y_i \\sim \\text{Bernoulli}(p_i)$（$n_i = 1$），当前模型（逻辑回归模型）预测概率为 $\\hat{p}_i = 1/(1 + e^{-\\hat{\\eta}_i})$，饱和模型中 $\\hat{p}_{\\text{sat}, i} = y_i$。\n",
    "- 当前模型对数似然：$\\ell = \\sum_{i=1}^n \\left[ y_i \\log(\\hat{p}_i) + (1 - y_i) \\log(1 - \\hat{p}_i) \\right]$\n",
    "- 饱和模型对数似然：$\\ell_{\\text{sat}} = \\sum_{i=1}^n \\left[ y_i \\log(y_i) + (1 - y_i) \\log(1 - y_i) \\right] = 0 \\quad (\\text{因 } y_i \\text{ 为 } 0 \\text{ 或 } 1)$\n",
    "- Residual Deviance：$\\text{Deviance} = -2\\ell = -2 \\sum_{i=1}^n \\left[ y_i \\log(\\hat{p}_i) + (1 - y_i) \\log(1 - \\hat{p}_i) \\right]$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcf72d1a-ea7e-4fa6-a46f-bfe4f8b2b238",
   "metadata": {},
   "source": [
    "AIC 的基本计算公式为：$$AIC = 2k - 2\\ln(L)$$其中，k 是模型中待估计的参数数量，L 是模型的最大似然值。在二分类问题中，我们通常使用逻辑回归模型来进行建模，下面以逻辑回归为例说明计算步骤。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c913d5bf-852a-4f87-8dc8-ec43eca2cc50",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = ynaffair ~ age + yearsmarried + religiousness + \n",
       "    rating, family = binomial(), data = Affairs)\n",
       "\n",
       "Coefficients:\n",
       "              Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)    1.93083    0.61032   3.164 0.001558 ** \n",
       "age           -0.03527    0.01736  -2.032 0.042127 *  \n",
       "yearsmarried   0.10062    0.02921   3.445 0.000571 ***\n",
       "religiousness -0.32902    0.08945  -3.678 0.000235 ***\n",
       "rating        -0.46136    0.08884  -5.193 2.06e-07 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 675.38  on 600  degrees of freedom\n",
       "Residual deviance: 615.36  on 596  degrees of freedom\n",
       "AIC: 625.36\n",
       "\n",
       "Number of Fisher Scoring iterations: 4\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating,\n",
    "                data=Affairs,family=binomial())\n",
    "summary(fit.reduced)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "41652a09-c9b9-46e6-8614-ac9b7d8b74c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analysis of Deviance Table\n",
      "\n",
      "Model 1: ynaffair ~ age + yearsmarried + religiousness + rating\n",
      "Model 2: ynaffair ~ gender + age + yearsmarried + children + religiousness + \n",
      "    education + occupation + rating\n",
      "  Resid. Df Resid. Dev Df Deviance Pr(>Chi)\n",
      "1       596     615.36                     \n",
      "2       592     609.51  4   5.8474   0.2108\n"
     ]
    }
   ],
   "source": [
    "print(anova(fit.reduced,fit.full,test=\"Chisq\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72e414ac-87dd-48d4-95a2-2d98ceda9dad",
   "metadata": {},
   "source": [
    "在进行 `anova(fit.reduced, fit.full, test = \"Chisq\")` 分析时，得到的结果是一个偏差分析表（Analysis of Deviance Table），用于比较两个嵌套模型（简化模型和全模型）的拟合优度。这里的 p 值是基于卡方检验（Chi - square test）计算得出的，下面为你详细解释其计算原理和步骤。\n",
    "\n",
    "**1. 卡方检验的基本原理**\n",
    "在比较两个嵌套的广义线性模型（如逻辑回归模型）时，我们可以通过计算两个模型的偏差（Deviance）之差来判断增加的自变量是否对模型有显著的贡献。偏差是衡量模型拟合优度的一个指标，偏差越小，模型拟合得越好。\n",
    "\n",
    "设简化模型的偏差为 $D_1$，自由度为 $df_1$；全模型的偏差为 $D_2$，自由度为 $df_2$。两个模型的偏差之差 $\\Delta D=D_1 - D_2$ 服从自由度为 $\\Delta df=df_1 - df_2$ 的卡方分布，即：\n",
    "$\\Delta D\\sim\\chi^2(\\Delta df)$\n",
    "\n",
    "**2. 计算步骤**\n",
    "\n",
    "**步骤 1：确定简化模型和全模型的偏差及自由度**\n",
    "从你给出的结果来看：\n",
    "- 简化模型（Model 1）：`ynaffair ~ age + yearsmarried + religiousness + rating`，残差自由度（Resid. Df）为 $df_1 = 596$，残差偏差（Resid. Dev）为 $D_1 = 615.36$。\n",
    "- 全模型（Model 2）：`ynaffair ~ gender + age + yearsmarried + children + religiousness + education + occupation + rating`，残差自由度为 $df_2 = 592$，残差偏差为 $D_2 = 609.51$。\n",
    "\n",
    "**步骤 2：计算偏差之差和自由度之差**\n",
    "- 偏差之差：$\\Delta D=D_1 - D_2=615.36 - 609.51 = 5.8474$\n",
    "- 自由度之差：$\\Delta df=df_1 - df_2=596 - 592 = 4$\n",
    "\n",
    "**步骤 3：计算 p 值**\n",
    "p 值是在原假设（增加的自变量对模型没有显著贡献，即简化模型和全模型的拟合效果没有显著差异）成立的情况下，得到比当前观察到的偏差之差 $\\Delta D$ 更极端值的概率。由于 $\\Delta D\\sim\\chi^2(\\Delta df)$，所以 p 值可以通过卡方分布的上尾概率来计算：\n",
    "$p = P(\\chi^2(\\Delta df)>\\Delta D)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6337f294-af04-415c-8831-d1cfac3073d3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 0.2108381\n"
     ]
    }
   ],
   "source": [
    "# 偏差之差\n",
    "delta_D <- 5.8474\n",
    "# 自由度之差\n",
    "delta_df <- 4\n",
    "\n",
    "# 计算 p 值\n",
    "p_value <- 1 - pchisq(delta_D, delta_df)\n",
    "print(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "386cde49-e1b0-4062-a0b7-e0712c4adad0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (Intercept)           age  yearsmarried religiousness        rating \n",
      "   1.93083017   -0.03527112    0.10062274   -0.32902386   -0.46136144 \n"
     ]
    }
   ],
   "source": [
    "print(coef(fit.reduced))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ee5a03bb-4376-46e5-b9de-ab40940470e7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (Intercept)           age  yearsmarried religiousness        rating \n",
      "    6.8952321     0.9653437     1.1058594     0.7196258     0.6304248 \n"
     ]
    }
   ],
   "source": [
    "print(exp(coef(fit.reduced)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f1889e39-a1ef-4842-9f2a-39102e4143a4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Waiting for profiling to be done...\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                    2.5 %       97.5 %\n",
      "(Intercept)    0.75404303  3.150622807\n",
      "age           -0.07006400 -0.001854759\n",
      "yearsmarried   0.04388142  0.158562400\n",
      "religiousness -0.50637196 -0.155156981\n",
      "rating        -0.63741235 -0.288566411\n"
     ]
    }
   ],
   "source": [
    "print(confint(fit.reduced))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a641264e-db70-4e81-b142-2fbfe523aaee",
   "metadata": {},
   "source": [
    "通过控制变量的方法，查看单一的特征发生变化是，预测概率的情况。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0639e454-db6c-4f36-9be8-96d1bd84bfba",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  rating      age yearsmarried religiousness\n",
      "1      1 32.48752     8.177696      3.116473\n",
      "2      2 32.48752     8.177696      3.116473\n",
      "3      3 32.48752     8.177696      3.116473\n",
      "4      4 32.48752     8.177696      3.116473\n",
      "5      5 32.48752     8.177696      3.116473\n"
     ]
    }
   ],
   "source": [
    "testdata <- data.frame(rating=c(1:5),age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried) , \n",
    "                       religiousness=mean(Affairs$religiousness))\n",
    "print(testdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b8bd23ef-d3a5-4d51-8470-80dd2ee3a928",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  rating      age yearsmarried religiousness      prob\n",
      "1      1 32.48752     8.177696      3.116473 0.5302296\n",
      "2      2 32.48752     8.177696      3.116473 0.4157377\n",
      "3      3 32.48752     8.177696      3.116473 0.3096712\n",
      "4      4 32.48752     8.177696      3.116473 0.2204547\n",
      "5      5 32.48752     8.177696      3.116473 0.1513079\n"
     ]
    }
   ],
   "source": [
    "testdata$prob <- predict(fit.reduced,newdata=testdata,type=\"response\")\n",
    "print(testdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2eb694f-41f1-4e42-ae51-7512e9f6c023",
   "metadata": {},
   "source": [
    "过度离势会导致奇异的标准误检验和不精确的显著性检验。通过检测过度离势增强模型使用的信心"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ceaf0d8e-1a17-4e19-a772-c70c9afc89ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.03247958314398"
      ],
      "text/latex": [
       "1.03247958314398"
      ],
      "text/markdown": [
       "1.03247958314398"
      ],
      "text/plain": [
       "[1] 1.03248"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "deviance(fit.reduced)/df.residual(fit.reduced)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e4950193-9a9a-4b55-9e92-9eba65c17adf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.340121988435398"
      ],
      "text/latex": [
       "0.340121988435398"
      ],
      "text/markdown": [
       "0.340121988435398"
      ],
      "text/plain": [
       "[1] 0.340122"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness + rating,\n",
    "                data=Affairs,family=quasibinomial())\n",
    "pchisq(summary(fit.od)$dispersion*fit.reduced$df.residual,fit.reduced$df.residual,lower=F)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faaa5b03-e2d9-40ab-b13b-a8e04387fd34",
   "metadata": {},
   "source": [
    "在原假设 $H_0:\\phi = 1$（不存在过度离势）成立的情况下，统计量 $\\phi\\times df_{residual}$ 服从自由度为 $df_{residual}$ 的卡方分布，其中 $df_{residual}$ 是标准逻辑回归模型的残差自由度。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "675e8817-676b-4111-b8c9-7ddc100c2370",
   "metadata": {},
   "source": [
    "## 泊松回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "2050eb36-f459-4c7f-b977-df6f5a615701",
   "metadata": {},
   "outputs": [],
   "source": [
    "# options(repos = c(CRAN = \"https://mirrors.aliyun.com/CRAN/\"))\n",
    "# install.packages(\"robust\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "a42dddb6-d17c-474d-a634-2d9d08fdc34c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data(breslow.dat,package=\"robust\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "06a2d9c0-4da2-4a33-8a74-546384f0e5df",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>'ID'</li><li>'Y1'</li><li>'Y2'</li><li>'Y3'</li><li>'Y4'</li><li>'Base'</li><li>'Age'</li><li>'Trt'</li><li>'Ysum'</li><li>'sumY'</li><li>'Age10'</li><li>'Base4'</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'ID'\n",
       "\\item 'Y1'\n",
       "\\item 'Y2'\n",
       "\\item 'Y3'\n",
       "\\item 'Y4'\n",
       "\\item 'Base'\n",
       "\\item 'Age'\n",
       "\\item 'Trt'\n",
       "\\item 'Ysum'\n",
       "\\item 'sumY'\n",
       "\\item 'Age10'\n",
       "\\item 'Base4'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'ID'\n",
       "2. 'Y1'\n",
       "3. 'Y2'\n",
       "4. 'Y3'\n",
       "5. 'Y4'\n",
       "6. 'Base'\n",
       "7. 'Age'\n",
       "8. 'Trt'\n",
       "9. 'Ysum'\n",
       "10. 'sumY'\n",
       "11. 'Age10'\n",
       "12. 'Base4'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] \"ID\"    \"Y1\"    \"Y2\"    \"Y3\"    \"Y4\"    \"Base\"  \"Age\"   \"Trt\"   \"Ysum\" \n",
       "[10] \"sumY\"  \"Age10\" \"Base4\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "names(breslow.dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "e8aefa7e-da81-40f6-b418-f327e3a51628",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   ID Y1 Y2 Y3 Y4 Base Age     Trt Ysum sumY Age10 Base4\n",
      "1 104  5  3  3  3   11  31 placebo   14   14   3.1  2.75\n",
      "2 106  3  5  3  3   11  30 placebo   14   14   3.0  2.75\n",
      "3 107  2  4  0  5    6  25 placebo   11   11   2.5  1.50\n",
      "4 114  4  4  1  4    8  36 placebo   13   13   3.6  2.00\n",
      "5 116  7 18  9 21   66  22 placebo   55   55   2.2 16.50\n",
      "6 118  5  2  8  7   27  29 placebo   22   22   2.9  6.75\n"
     ]
    }
   ],
   "source": [
    "print(head(breslow.dat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "6b048a97-0046-4768-a456-cf2987561ec5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      Base             Age               Trt          sumY       \n",
       " Min.   :  6.00   Min.   :18.00   placebo  :28   Min.   :  0.00  \n",
       " 1st Qu.: 12.00   1st Qu.:23.00   progabide:31   1st Qu.: 11.50  \n",
       " Median : 22.00   Median :28.00                  Median : 16.00  \n",
       " Mean   : 31.22   Mean   :28.34                  Mean   : 33.05  \n",
       " 3rd Qu.: 41.00   3rd Qu.:32.00                  3rd Qu.: 36.00  \n",
       " Max.   :151.00   Max.   :42.00                  Max.   :302.00  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "summary(breslow.dat[c(6,7,8,10)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74cfa171-e4dd-4c4e-a84b-b721cb8967aa",
   "metadata": {},
   "source": [
    "- `sumY`：随后八周内癫痫发病数，作为响应变量；\n",
    "- `trt`:治疗条件；\n",
    "- `age`:年龄；\n",
    "- `base`:前八周内的基础癫痫发病数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "014785a8-4427-4783-bf64-19ca2a88ba89",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        Trt Age\n",
      "1   placebo  28\n",
      "2 progabide  31\n"
     ]
    }
   ],
   "source": [
    "\n",
    "library(plyr)\n",
    "print(aggregate(breslow.dat[c(\"Age\")],by=list(Trt= breslow.dat$Trt),FUN=length))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c5ac111-ad09-43d6-a9cc-5736b5576902",
   "metadata": {},
   "source": [
    "用length进行计数……placebo 就是安慰剂。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "5a1877f7-3f46-440e-906d-7694b1f9bf89",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEX9/v0AAABMTUxnaGd7\ne3uLjIuZmpmmpqaxsrG7vLvFxsXOz87T09PX2Nff4N/n6Ofu7+79/v03sFGxAAAAEnRSTlP/\n/////////////////////wDiv78SAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3d\n2YKiyhJA0ZMMosVV4f9/9sqgZKAiBGGa0Hs9dNvVSFGku5RB+a8GsNp/v14AYA8ICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQZChOR6WXF9fOHFZIene72fdnyvT1NNOSbO+d+7PNy+kBaX\nF5Ou+TZbdC2y24+c5H9G8zsfUsvZ3cUwLCFDuinvX3ia6Jy4p3u9m/b5XitW5bFZLi+k4r6s\nx+dpYxixgI6PcUtKi/nlj9mdLWb3EMOwBA7JnfsvvJro5RemV9L9f1esytvvSP/J5zRa1umF\n3LXsaeDs5vfq+V4thmEJFFLzV1U2z+sfJpr99Tn/O8toFrdlLKrbi5rbqOdrZ71tzfPRoenn\n0txKV8/v9nyUnJpV2/yu2t26DRhS3T5Ky+EL1bH5HdW+Yu5/T3W3rqkrxDPSKXXJ4Srm1f4t\n7tV9h/LQDH85fOPy9j0Oo99/TxO9WNjqcaNIXFL43334vTpeIG/xxf28n3Qzro9X4rcVlh7b\nH2T46V6vQzlmrh+zzuXWUf/vSz+755nc7pXe2j0lLjvXr2b117w8TAu5NM+Pp5cz9x4K9gMS\nNqSy2xrpvnBN+odjNkoivX+pv3P24lXhm5Durx7y+0TFi1cS/kTDLHq3b577L2PuS+l99w8h\ndYvv38//STfj9jRUjL/2+OmeV7T/d/P19n8Tr6TbSJxGs3uaSfeFayHWuD+rx6vDc/30UBmt\n5amHwhcGJGxIVfcSofvC7RfG7TdClTUrWDwwXfsfw6D0En9e7mVIj63ZbvUN9/X3yomJhln0\n2m2k5PB3fwjcV7n33YfZjhfIW3z/fv5PuhnZiw2Zx0/3vKJr7+9hBaViflc5t3ejlbwcwmZW\nt8HJqrbJ4Zfg36vH04eHwhcGJGxIwyOx+7NZsyKu7vdS5d3LdbuMyqT92V+9iBj+vj3hudur\n8Kp5TV8O9z2IUl5MJJb2PgJp+6KgG7rqIL97o5B3HxaoXXxxP/mTbsR4vfRfa3+6d+vQG7Nz\nN2bl+/m9mMlt5s0vsvTS/vU8q7SP8fEA8h8qYi1PPxS+MCA/DKn5xXMoxxO50T7y+7/Lx2+h\n0f96/zzcf8MUj5eQzX0rMYYvJpKL2+wTaTWvbG5Z9WOVy2n7jl4sUPtlcT/5k27E42ft10Z/\ns/0x3q3DYbJ266P0Xwo8recXMzmLv97OangY+Q8VsZanHwpfGJAfhtQfpTjIY0uufwSOUxF3\nfRPS477XlxM/5vNyIs/179C+vDiNX1kO0xbjY2LDAt1/SQ73kz/pRrwJqfL/frOivbumz/Pz\nvvBqJi//GmZ1/WsPEovhfvF4ml7CLwzID7eRHoc+2+3Ipwe9LqR6YuLx0rwNqXHN22UdghDT\nnh5b4u8XyLuf+Ek3Ir8fO5Ih1f7fH0ZFrNt0vM31ZiZvQ2pu/KXeWh1/G38tf1hC+wH54V67\nW1d/3a6VrH4xRqNfgkufkZJX86wnJuok9//t/yMR/3ufthz2aL1ZIHk/8ZNuxPHF6zLvh51a\n0aMx6zzttXszk6e/vFn93f5MD6fL65D8tfzpoWA+IGFDyp6378uDe1WIvybvL5L7rZRmBZ1f\nhpRPv3TvvJuoc3g8fLrXBLm/vXyf1uvo3QLJ+4mfdCOazYnxFqz3043X4WgltP9d+gdez8Pv\n/3N7HOnNQDz95c0qFRumL5/4+rU856FgOiAhQzrnYi9y+njBff+VUdUvQ2pKanbbNCe/Je1j\n+Jx4q3u414edSZ3pvXbN/7YH/8ruVNa/bq/Rn/+06Xf0aoEa4n7yJ92KZgXl3ZkN459uvA6f\nV8J9zLwnoceZDe2M3w7E01/erPovvn5GEmt5+qHwhQEJfa6d9wvl9rNm18eBv2YQivp1SJ2k\n+dkPw7/r53sNZ3MdvFmMWnkzUW84/tDtIr0fD/J2JYlFGC3QY27+/eRPuhniXDtxvGi8Dp9X\nwn3MvNlViff16/NM3LuQhlll7RosE/cqJLmWJx8KXxiQX579fd/ia1+nHpzcVvLWZO7d99rd\nLvr/Hd/rvvrkmyucbOX1ROP/7c94Lvt/FcO04lEzWqDH3MT9xE+6HcPZ3+n4tH25Dp9XQvcT\ny2356rFq0+uLmbg3IXmzOt+Hpv39NH6oyLU8+VCwH5Cw70eqHl9o/25fpGb9s38+2mTxVugp\nG062utzuk/35r9bFvdo3E705Bezh5UTyf11+7Be2KpqzhvwHkghptEDD3Pz7yZ90O6pj3o7b\n/Zwpb12JdfhiJfylLrmPdy3ulA57nV8NxPNf3qya75McLtf7WSlyseRannwomA8I75CFuaff\nTVHM6rsICeYICTBASIABQgIMEBIAFUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBAS\nYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBAS\nYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBAS6vYKk/AsX4OEhJuNXKg1FEKCDiEJhAQdQhII\nCTqEJBASdAhJICToEJJASNAhJIGQoENIAiFBh5AEQoIOIQmEBB1CEggJwvmYt6eO5cV5ekJC\nEggJnir1TsPMJiclJIGQ4Clc8ndpb13LxBVTkxKSQEjwJO7yuH1xydSkhCQQEjzibTXT77Eh\nJIGQ4OEZSYuQ4LltI5XX9hbbSMsQEnyZt9curaam3EVIqreIv57T8rsQ0p6di/Y4UpIf938c\nqa3IKCVCgs4eQvL+tJnXIoSEeg8hudHfJjNbIGhIBp96hCX+oVOE/q2Q/idtf/ii9k+dIkRI\n+JZ/6xShf2kbiZCC+rcOyP5Le+0IKagPpwjtbmv13zmOREhB/VvPSJYICR5OEdIiJPj+sVOE\n7BAShH/pFCFLhAQdRkIgJOgwEgIh4YVT4tLT9CSMhEBI8F1yl5zq4z9xipApQoLn0hZUuENV\nX3M3+ZzESAiEBM+hOXZUdEdiK5dOTcpICIQET3fGjMu9f7ydNMDibAghwdO189e9puMUoQUI\nCZ5Ds3XUqQ6cIrQAIcFTJY/Xc276CYmQJEKCUNzzSSafjwhphJCgw0gIhAQdRkIgJOgwEgIh\nQYeREAgJOoyEQEjQYSQEQoIOIyEQEnQYCYGQoMNICIQEHUZCICToMBICIUGHkRAICTqMhEBI\n0GEkBEKCDiMhEBJ0GAmBkKDDSAiEBB1GQiAk6DASAiFBh5EQCAk6jIRASNBhJISgIZ2P7dXg\nXF58uBrcgJBixUgIAUOqUu/6pNOXDBkQUqwYCSFgSIVL/rqLz3+8YvaAkGLFSAgBQ0rc5XH7\n8uHzcB8IKVaMhBAwJHGVkOlLhnjTEVKkGAmBZyToMBJC2G2k8treYhtpBxgJIeTu78zba5dW\nn6dvEFKsGAkh7HGkoj2OlORHjiNtHiMhcGYDdBgJIZ6QnG/4KiFFipEQfhDSKXHpaXoSQoof\nIyGEDOmSu+RUH2ecIkRI8WMkhIAhXdqCiua62dfcTT4nEVL8GAkhYEjt5eaL7khs5dKpSQkp\nfoyEEPwUIZd7/3g76XCLkCLFSAjBQ/rrXtNNnyJESPFjJISgL+0O99MZqsP0KUKEFD9GQgj5\nxr7k8XrOfThnlZDix0gIQY8jFfd8kg+nrBJS/BgJIZ4zG3yEFD9GQiAk6DASAiFBh5EQCAk6\njIRASNBhJARCgg4jIRASdBgJgZCgw0gIhAQdRkIgJOgwEgIhQYeREAgJOoyEQEjQYSQEQoIO\nIyEQEnQYCYGQoMNICIQEHUZCICToMBICIUGHkRAICTqMhEBI0GEkBEKCDiMhEBKE87G9qqLL\niw9XVWQkBEKCp0q9q73NvfQOakKCULjk79Le+njleUZCICR4End53L7MvdABakKCIK62M/fS\nO6gJCQLPSFqEBM9tG6m8trfYRlqGkODLvL12aTU1JSMhEBKEc9EeR0ryI8eRliAk6DASAiFB\nh5EQCAnPTqlzeTk9DSMhEBI83aGjfo/D9PVJGQmBkOBpQypcUdX1tXCnyUkDLdJGEBI8bUiJ\na/d7Vy6dnDTIAm0GIcHThnQ/Nej5FCHnC7xokSMkeNo8DveQOEVoPkKCx7n8eCrd3+1mVXCK\n0AKEBI/3ss25hFOE5iMk+C6X0ynP210OxWRHhCQREnQYCYGQoMNICIQEX3VwLutPDuIdsgsQ\nEjxV0n0WV/sPQlqAkOBpTwuqTkn7SVyEtAAhwZN06/eapFdCWoSQ4Lm3U2UZIS1CSPCk7n7w\nKM0IaQlCgufkDv2tq8sIaQFCgq941FN+OMGbkRAICcIlv9+6HghpPkKCDiMhEBJ0GAmBkKDD\nSAiEBB1GQiAk6DASAiFBh5EQCAk6jIRASNBhJARCgg4jIRASdBgJgZCgw0gIhAQdRkIgJOgw\nEgIhQYeREMKHtOxqcIQUK0ZCCBiS6mpwhBQrRkIIHdLSq8ERUqwYCSF0SEuvBkdIsWIkhNAh\nvb8anJh0uEVIkWIkhNAhLb0aHCHFipEQgoakuBocIcWKkRCChqS4GhwhxYqREEIeR9JcDY6Q\nYsVICJzZAB1GQognJOcbvkpIkWIkhJAhaa4GR0ixYiSEgCGprgZHSLFiJISAIamuBkdIsWIk\nhIAhqa4GR0ixYiSE4Gd/L7waHCHFipEQAoakuhocIcWKkRAChqS6GhwhxYqREELu/tZcDY6Q\nYsVICEEPyCquBkdIsWIkhHjObPARUvwYCYGQoMNICIQEHUZCICToMBICIUGHkRAICTqMhEBI\n0GEkBEKCDiMhEBJ0GAmBkKDDSAiEBB1GQiAk6DASAiFBh5EQCAk6jIRASNBhJARCgg4jIRAS\ndBgJgZCgw0gIhAQdRkIgJOgwEgIhQYeREAgJOoyEQEjQYSQEQoIOIyEQEnQYCYGQoMNICIQE\nHUZCICToMBICIUGHkRAICcL5mHfXzC7O0xMyEgIhwVOlbpBNTspICIQET+GSv0t761omrpia\nlJEQCAmexF0ety8umZqUkRAICR5xRdK5lydFTUgQeEbSIiR4bttI5bW9xTbSMoQEX+bttUur\nqSkZCYGQIJyL9jhSkh85jrQEIUGHkRAICTqMhEBIEDhFSIeQ4OEUIS1CgodThLQICR4OyGoR\nEjwfThFyvoCLtQGEBA/PSFqEBA+nCGkREnycIqRESBA4RUiHkKDDSAiEBB1GQiAk+KqDc1nZ\n3eYdsgsQEjxV0p1o1/6DkBYgJHgKd7rVdEra0+wIaQFCgifp1u81Sa+EtAghwXNvp8oyQlqE\nkOBJ3f0gbJoR0hKEBM/JHfpbV5cR0gKEBF/xqKf8cII3IyEQEoRLfr91PRDSfIQEHUZCICTo\nMBICIUGHkRAICTqMhEBI0GEkBEKCDiMhEBJ0GAkhaEiKj8MlpFgxEkLAkFQfh0tIsWIkhIAh\nqT4Ol5BixUgIAUNSffggIcWKkRAChqS6YjYhxYqREHhGgg4jIYTdRlr+cbiEFCtGQgi5+1vz\ncbiEFCtGQgh7HGn5x+ESUqwYCYEzG6DDSAjxhPT6IlaEFCtGQggZkubjcAkpVoyEEPIUIc3H\n4RJSrBgJIejub8XH4RJSrBgJIegB2favZR+HS0ixYiSEH5witOjjcAkpVoyEEDAk1cfhElKs\nGAkhYEiqj8MlpFgxEkLI3d+aj8MlpFgxEkLQA7KKj8MlpFgxEkI8Zzb4CCl+jIRASNBhJARC\ngg4jIawMKT1ezRbFQ0jxYySElSE1b9H7QkuEFD9GQlgZUvV3+EZLhBQ/RkIw2EY6H1Prlggp\nfoyEYLOz4dK8Q+K0fmnuCCl+jIRgElKZzfgY4iUIKX6MhLA+pOp4ezpKy+pWU/76DssRUvwY\nCWFtSOdmZ0PRffLj9PlzSxBS/BgJYe1xpNuT0en+7ojpT09dgpDix0gIa48j5aXZovizHW4R\nUqQYCWHtcSSzBREIKX6MhLB2G6kqmtdzSWFbFCHFj5EQVoZ0Tdo9DM4lpuc2EFL8GAlhZUiZ\nOzTPRVVht+u7QUjxYySE1Setjm+YIKT4MRLCypCS/pOBKkL61zASwsqQCpc1V2g5Z9MXDluK\nkOLHSAhr99rdLx5md55dg5Dix0gIq8+1+2uuHZYZnvndIKT4MRICn9kAHUZCICToMBICIe3a\n4fMkWoyEsDak5m3mo6tVGiAkIy69fG3W35rxNq0M6fjisq8GCMlI7tzxS7NmJITVB2SN99d1\nCMnKybnsKx89SEiS1SlCtgjJzDUz/ViaASMhrAwpd195RxIhGbq9/M6+MEqMhLD6bRTtKULW\nCMlSVXx3Oxa1xUcWs7MhekdC+jpC2j1e2oXAAdmdY2dDGIS0b+z+DmR1SGXevKrLbQeLkIxw\nQDYUk/cjNZ8NyYefxIhThEJZGdLpth3bhHRypqdHEpKR5aNyPubtvqO8+HBcg5EQDD6zof9A\nLqslahDSj1Sptx92+l3PjIRgcIoQIe1H4ZK/7sXgtUymP4eDkRBWhpT2z0gXl5otUk1Itopk\n9rG+xA3bVJfpiyIwEoLNNlJpfBY4IRkqFhw0F5NMT89ICGv32uV8ilDslhyQ5RlJy+Q4ksv/\njBanR0iGlmy+3raRyu5ABttIy3Bmw+4VS97qknl77dLJ+zESAiHtX7bkJKFz0b5aT/Ijx5GW\nIKT9K798hj5q3kbxD/j6B9SgJqR/gOLQxCl1H68OzEgINi/tzpnpdcYIydKS33HdtP0eh+kL\njDASgtE2UsVJq9E6Lthr14ZUuOaSwNdi+pmMkRCsdjbw0i5ex/kfUNMOY3/1uGr6tC9GQjAK\n6TR9FHwpQjK0ZDtWnID8PL1zX9oo3j6znQ2mb8QkJEOLQzrcQ+IUofmMQkptP1+DkH7Eufx4\nKl1zxldVcIrQAhyQhcd76nIu4RSh+QgJvsvldMrzdpdDMb2zj5EQ7A7IWm5/EpKhrx80R01I\n/wBCCmHtS7tj0pxJck54Y1/svnf2CerVIR37d1RenOkoEdIXfO3sE9QmnyIkb5ggpG/gpd0X\nrQwpeTwj8SlCsZtz9smCbV6rkdjJORIrQ2re43/7i08RitiSs09OoUP6woci/sbanQ339/hP\nn3K/FCEZWnT2yWX2biOjkAzn9VOrD8j+tZ8i9OFNYEsR0s9c5v5ONBkJN/p7uzizAdLJzbuA\nBSEJ4UNa9jZmQlrvlNb1NXWp7WWzCUlYHdL8C42p3sZMSKuVzYpvP/7btCS2kQSTnQ31rAuN\nqd7GTEirZe6vPTzxZ/u50uy1E1aGtORCY6q3MRPSat3VQopYD8hyHKmx5EJjH97G/HqpCGm1\n9qW3K2MNaScMThFaFNLStzET0mqZu5TN6o7ypd1urAxpyYXGVG9jJqTVyv6kBudMj/YxEoLN\nNtKsU4RUb2MmpPVO3QVaUttr7zASwtq9dksuNKZ5GzMhxYqREEyOI3GhsX8QIyHEc4rQ69OO\nCSlWjISwMqTc9qzvO0KKHyMhWL1D1hYhxY+REAx2f38BIcWPkRBWhlTlS650oHj3JSHFipEQ\nVr+0m/+Zaaq3MRNSrBgJIWBIqrcxE1KsGAkh6O5vxduYCSlWjIQQ9jjS8rcxE1KsGAlhRUhf\nfB8JIcWPkRBWh/SVnAgpfoyEQEjQYSQEQoIOIyEQEnQYCYGQoMNICIQEHUZCWBXS7FN+9EtF\nSLFiJARCgg4jIcTzDlkfIcWPkRAICTqMhEBI0GEkBEKCDiMhEBJ0GAmBkKBjNRJcjeKLCCl+\nNiPB9ZG+ipDiZxSS4bx+ipCgYzISbvT3dhESdAhJICToEJJASNBhG0kgJOiw104gJOhwHEkg\nJOgwEgIhQYeREAgJOoyEQEjQYSQEQoIOIyEQEnQYCYGQoMNICIQEHUZCICToMBICIUGHkRAI\nCTqMhEBI0GEkBEKCDiMhEBJ0GAmBkKDDSAiEBB1GQiAk6DASAiFBh5EQCAk6jIRASNBhJITN\nhfTFC25iCVa9sLmQeIKKBKteICTosOoFQoJwPubta+a8OE9PyKoXCAmeKvW2P7PJSVn1AiHB\nU7jk79LeupaJK6YmZdULhARP4i6P2xeXTE3KqhcICR5xPGH64AKrXiAkeHhG0iIkeG7bSOW1\nvcU20jKEBF/m7bVLq6kpWfUCIUE4F+1xpCQ/chxpCUKCDqteICTosOoFQoLAKUI6hAQPpwhp\nERI8nCKkRUjwcEBWi5Dg+XCKEG9OfouQ4OEZSYuQ4OEUIS1Cgo9ThJQICQKnCOkQEnRY9ULQ\nkBRHzQkpVqx6IWBIqqPmhBQrVr0QMCTVUXNCihWrXggYkuoYBSHFilUvBAxJ9cEahBTUgk9W\nX7rqx5/aPsOKHyQ4npHgOX0vpO/P6KfCbiMtP2pOSGFdkundQANCEkLu/tYcNSekwC7Tv+IG\nhCSEPY60/Kg5IYV28l6BTyEkgTMboMOqF+IJ6fVWLiHFilUvxBOSj5Dix6oXCAk6bCMJhAQd\nQhKCntmgONhHSLEiJCFgSKqj5oQUK0ISQr600xw1J6RYEZIQdBtJcdSckGLFqhfC7mxYftSc\nkGLFqhfYawcdVr1ASNBhG0kgJOgQkkBI0CEkgZCgQ0gCIUGHkARCgg6rXiAk6LDqBUKCDqte\nICTosI0kEBJ0CEkgJOgQkkBI0CEkgZCgQ0gCIUGHVS8QEnRY9QIhQYdVLxASdNhGEggJOoQk\nEBJ0CEkgJOgQkkBI0NlFSHYXqiUk6Oxg1bcVGaVESNDZwap33p8281qEkFDvISQ3+ttkZgsQ\nEuo9bCMREiFFgJBezmwBQkK9h5DYRiKkCOwhpM9X6po/r+V3ISTUhDSe1/K7EBLqPey146Ud\nIUVg+6uenQ2EFIHtr3pCIqQIbH8biZAIKQLbD4ltJEKKwB5C4qRVQvq5HYTE2ygI6fd2EZId\nQoLOLlY9z0iE9Gs7WPVsIxHS7+1g1bPXjpB+b/vbSBxHIqQIENLLmS1ASKgJ6c3MFiAk1HsI\niW0kQorAHkJirx0h/dwuVj3HkQjp11j1AiFBh1UvEBJ0drCNZImQoENIAiFBh5AEQoIOIQmE\nBB1CEggJOqx6gZCgw6oXCAk6rHqBkKDDNpJASNAhJIGQoENIAiFBOB/z9loneXGenpCQhK2H\nNPb9Zdu1KvXWZTY5KSEJWw9p/N/fX7ZdK1zyd2lvXcvEFVOTsqoFQoIncZfH7YtLpiZlVQuE\nBI94bTz9QplVLRASPL94RtrHmBESPLdtpPLa3gq3jbSPMSMk+DJvr11aTU1JSAIhQTgX7XGk\nJD/u6TjS02GSz5Z+h+ULRUioNxbS9785IWG+/R37JiRC+qJT6lxeTk+zj1VNSIT0Dd3zTL/H\nYXKnHSGtnhEh7VcbUuGKqq6vhTtNTmr2Pa1m9FPhQ1r2soGQgmpDSly737ty6eSkZt/TakY/\nFTAk1csGQgpKXJ9Bc4qQYj/z93dNhxA6pKUvGwgpqHaMDveQFKcIjQfkKwxHeYvbSKqXDYQU\n1O1V9/FUur/bzarQnCJESPOtCun9y4bXz92EFJS3/p1LFKcIEdJ8q0Ja+rKBkMK6XE6nPG9f\nOxSTHRHS6hnpQ1K8bCCkWBHS2hnpQ1K8bCCkWO0jJDMhjyNpXjYQUqwISeDMBugQkkBI0NlH\nSFvcRlqAkOJHSGtnREioCWn9jAgJNSGtnxEhoSak9TMiJNR7CckMIUGHkARCgg4hCYQEnX2E\nxDbSy38SUjiEtHZGhISakNbPiJBQE9L6GUUd0gY+JmMnCGntjKIO6XurHNI+QjJDSNAhJIGQ\noENIAiFBZx8hsY308p+EFA4hrZ0RIaEmpPUzIiTUhLR+RoSEmpDWz4iQUO8lJDOEBB1CEggJ\nOjGG9M1rMn0464yQoBNlSL/75oQEHUKasTomERJqQpq1OiYREmpCmrU6JhESakKatTomERJq\nQpq1OiYREmpCmrU6JhESakKatTomERJqQpq1OiYREmpCmrU6JhESakKatTomERJqQpq1OiYR\nEmpCmrU6JhESakKatTomERJqQpq1OiYREmpCmrU6JhESakKatTomERJqQpq1OiYREmpCmrU6\nJhESakKatTomERJqQpq1OiYREmpCmrU6JhESakKatTomERJqQpq1OiYREmpCmrU6JhESakKa\ntTomERLqtyH98FODCekJIcWPkGasjkmEhDrSkH73zQkJOlFuIxGSREjxI6QZq2MSIaGONKTf\nfXNCgk6U20jG0bxCSDAV46olpBFCip/ZqjUcI0IaIaT4EdLan4KQUBPS+p+CkFAT0vqfgpBQ\nRxrS73YZEhJ09rFqf/nrYFMhTf820f3y+dYP/dPvvljkizcTIc0MafoJauHkX/6hxwsT9Lsv\nFvnizURIhPRjMW4j/RAhEZIOIQmEREg6hCQQEiHp7CMktpEI6ceWLp7ZERtThERIPxb54s1E\nSIT0Y5Ev3kyEREg/FvnizURIhPRjkS9eaIRESDqRL15ohERIOpEvXmhBQzof83a3Zl6cpyck\npJ+FpBijLdviNlKVeocIsplLRUhBqcZoy7YYUuGSv0t761omrpialJB+FJJqjLZsiyEl7vK4\nfXHJ1KSE9KOQVGO0ZVsMSZzy8Xz+x+uTQxQnluyGcj2voRqjDfj+2o/zGQk/8s89I5kJu41U\nXttbH19/40cWjBEhCSF3f2feE2daaeeCb5o/RoQkhD2OVLTHKJL8+OEYBX5m9hgRkhDnmQ2I\nHyEJhAQdQhIICTqEJBASdAhJICToEJJASNAhJIGQoENIAiFBh5AEQoIOIQmEBB1CEggJOoQk\nEBJ0CEkgJOgQkkBI0CEkgZCgQ0gCIUGHkARCgg4hCYQEHaP5dskAAAecSURBVEISYg9J8TlK\nmxVyva7367UVm+VrMGxIv5jsN991YyGZ+enP/ctvTkg/nWx3CCkEQto9QgqBkHaPkEIgpN0j\npBAIafcIKQRC2j1CCoGQdo+QQiCk3SOkEAhp9wgpBELaPUIKgZB2j5AA6BESYICQAAOEBBgg\nJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAgXAhFYlLimpiAv/zyz9O\n/MHJ2c3zlN7vZ7mEWzfng+Z/dSkB7/u+vvkFwULK2sdg+n6Ci/cw/TjxB5d+pVnMs2jvl1S2\nS7h5hCSFCunskkt9Sdz57RQXl8+feNrtrveQVs/z4g5V8wx3MF3C7dtISB+/aiVUSIUrb3/+\nuePbKU7D/32eeNLJZf1KM5hn3s2pmaHdEu4AIUmhQsrdtRa/05+c3Gn+xJNcUT9CMpuns5zb\nVt3WQnHbNOxv1nWZu/6f7UZjdm1vNRuVp9Hk7VfT0/Msl33PKm1XuTev27ctxgszXkxvomHh\njIUKyTn/r1dyVx7ET6/+FXJ53NVsnpXLLJdwq5w7NpuGWd396Mdum7FdI9l9S/K2moZp8vvN\nfqMyW/k98/a7efNqbx5GC+N933aE2ony9mauXYxPYgrJiZ9+1cP0EZLRPE/NSznTJdwk128a\n/nU/umtu/LUr4c9lVX1oHsZlc6vKmhXmTf433FzzPbOqFvMq+5tyYUaL+TdM5C2csXhCaldE\nVTQvn8xCsprnNcmtl3CTXPsILJuXV8OP3v+mPzfP20lzq3mwV900j8nz/uby5wLxPdudO968\n7jflwowWs122biJv4YzFE1KnanYpm4VkNM8q8cbfaAk3yfuxu5vX8pj1j91hmsdRgqfJNSvs\n7UyevjJemPsE4p/qqy1/ECqkZO6KbKaYPfHkXOzmmYnjRTZLuEnjh3D2eFRGENLTwuwypG43\n1/Xzc2r3BDxz4sm5WM3zmvZ7o0yXcJNGj9CDS0/l9Smkd5PXXw3peWFehrT4+88UKqRj/7K1\neDtF0r58bR+bnyf+pF9hFvMcXtebLuEmuX5r4+A9LK/908GwjVSOJ/e2kZb/5hl9z1rMy9tG\n8hbG+77db71movNw8xviObOhaB6VVXukc/15A/0qN5jnddg+Nl3CTbrvDivvIZ3rS7dZcmp2\nh7UrqN1Hdvt3LiZfv9eufIzqu712j4UZ3cWbyFs4Y8HOtUs/7cCvkuGgxMeJP+lXucE8D8Or\natMl3CTnukMydbeGi37VNL9OhuNI3bZKcm2mOdwnX3McSXzP0byyx+h4C+N9336X4v1gk7dw\nxoKFVLWnS3+aoj9e/XniD+6rfP08/c1TyyXcpNtayPs10K6Q24MzO/ev124P5PxxZoM7XLtp\nCpf0p1GdEvWZDeJ7jubVnFBxHi2M9327uxzFmQ3dwhnj/UiY7xd7Ked8zwheGhAS5osupO4Y\neR7BTh9CwnzRhdSfYJcEW5q3CAnzRRdSfcqcS3//fERIgAlCAgwQEmCAkAADhAQYICTAACEB\nBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEB\nBggJMEBIgAFCAgwQEmCAkAADhAQY2FpIVZE6l42v/Lb8KgmXQ+IO8y7M+63L96LhPI8vbnCV\nbyyk/jKu3cVKB4tD6i84ms64BmL6g0uZ/ENehbTFVb6xkA4uuz32r9nKS7QdXXL7pVcd51yV\n9xfXBPrXjNbxFlf5xkJyrn0qqtat6us9oIM7zPiWa74V5iCk0OQqPqUueVzu+vH6oJum+1KV\ntpfcfkzYKVx/pe0qP/Xz8S+b3d31mrcXxpYv3vEV9zXcDdgmV/nGQir8S7t3a7y5oPXbkG6T\nFP6Encxd/JlmYj73u7ZbY0dCCmEIqRmwTa7yjYVUt5cMPbc3S5dVdZW50nuium1DiRqySk7Y\nkcP055JLfUma62OP7npy6TZfZ2zNEFI7YFtc5VsLqS4PzU67pom83V6qmtdu9zWfNc8rfg1t\nct6EHTlQeVtY+equ3a0AP9U/bgjp7P9zSzYX0s35mDQr3Ntr2q/57On1Wd3dHO1elQP1mOr5\nroQUxhCS+OeWbDGkur50r7lGIWXdPrjPIeWPbaSyIqQIEFJgjzU8eoQ3N69Jf3DpRUij2Rzv\ne+3O3kYQIf0OIQWWu24nduWS+7ZN67bmHx11w3D2GvAm7DyOI2XN/O7bSPnzXQkpDEIK7PYY\nP91ejJ3bANq9bfWpD+DRUZ3e/rPKvAa8CXuHdndFc6SoFnvtxnftjyiF/jH/OU8hbW+Vbyuk\n+zly3TGh7vhP8+TSHj26bwidmj9z/8lkmPAu88+1G44jje/a/Jm6Njd80SikLa7yjYXUnrXt\nsr/uH6fbGm8P0IqQ6mPiDvJV2WPCh798mE19SvozG8Z3bf48p9sb1a0ZhbTFVb61kIAoERJg\ngJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJg\ngJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTDwfxYwtw1iOwiZ\nAAAAAElFTkSuQmCC",
      "text/plain": [
       "Plot with title \"Group Comparisons\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# opar <- par(no.readonly=TRUE)\n",
    "par(mfrow=c(1,2))\n",
    "with(breslow.dat,{\n",
    "    hist(sumY,breaks=20,xlab=\"Seizure Count\",main=\"Distribution of Seizures\")\n",
    "    boxplot(sumY~Trt,xlabel=\"Treatment\",main=\"Group Comparisons\")\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e2dce63-fe2f-487b-8cd5-3c67eadf5b00",
   "metadata": {},
   "source": [
    "within 函数用于在数据对象内部进行修改，它会返回一个修改后的数据对象，而原始数据对象不会被改变。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "914e4c43-7180-4e3f-a818-a60b2e92648d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   name salary bonus total_income\n",
      "1 David   5000   500         5500\n",
      "2   Eve   6000   600         6600\n",
      "3 Frank   5500   550         6050\n"
     ]
    }
   ],
   "source": [
    "# 创建一个数据框\n",
    "employee_data <- data.frame(\n",
    "  name = c(\"David\", \"Eve\", \"Frank\"),\n",
    "  salary = c(5000, 6000, 5500),\n",
    "  bonus = c(500, 600, 550)\n",
    ")\n",
    "\n",
    "# 使用 within 函数计算总收入（工资 + 奖金）并添加到数据框中\n",
    "new_employee_data <- within(employee_data, {\n",
    "  total_income <- salary + bonus\n",
    "})\n",
    "\n",
    "print(new_employee_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "7733069b-6d34-4878-bace-ff73ec2d85fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = breslow.dat)\n",
       "\n",
       "Coefficients:\n",
       "               Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)   1.9488259  0.1356191  14.370  < 2e-16 ***\n",
       "Base          0.0226517  0.0005093  44.476  < 2e-16 ***\n",
       "Age           0.0227401  0.0040240   5.651 1.59e-08 ***\n",
       "Trtprogabide -0.1527009  0.0478051  -3.194   0.0014 ** \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for poisson family taken to be 1)\n",
       "\n",
       "    Null deviance: 2122.73  on 58  degrees of freedom\n",
       "Residual deviance:  559.44  on 55  degrees of freedom\n",
       "AIC: 850.71\n",
       "\n",
       "Number of Fisher Scoring iterations: 5\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit <- glm(sumY~Base+Age+Trt,data=breslow.dat,family=poisson())\n",
    "summary(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "a9bde572-81ff-4db6-8af1-095619d35f07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " (Intercept)         Base          Age Trtprogabide \n",
      "  1.94882593   0.02265174   0.02274013  -0.15270095 \n"
     ]
    }
   ],
   "source": [
    "print(coef(fit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "48e507fe-17b3-41c3-967a-e38c45112538",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " (Intercept)         Base          Age Trtprogabide \n",
      "   7.0204403    1.0229102    1.0230007    0.8583864 \n"
     ]
    }
   ],
   "source": [
    "print(exp(coef(fit)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9debbe81-bd3e-4e98-8d10-73e9ad5b428f",
   "metadata": {},
   "source": [
    "过度离势（Overdispersion）是指数据的离散程度比预期的理论模型更大的情况。在以下几种情况下，通常需要计算过度离势：\n",
    "- **数据分布与理论模型不符**：当观测数据的方差明显大于基于某种理论分布（如泊松分布、二项分布等）所预期的方差时，可能存在过度离势。例如，在对某地区疾病发生次数进行建模时，若使用泊松分布假设，理论上均值和方差应相等，但实际数据中方差远大于均值，此时需考虑过度离势，可能是因为存在未考虑到的因素导致疾病发生的个体差异或环境因素等影响。\n",
    "- **模型拟合效果不佳**：如果使用常规的统计模型进行拟合，得到的模型对数据的拟合效果较差，残差分布呈现出非随机的模式，如残差的方差随预测值的增大而增大，这可能暗示存在过度离势。此时计算过度离势有助于评估模型的适用性，并考虑是否需要对模型进行改进，如采用更复杂的模型或对数据进行变换。\n",
    "- **分析数据的异质性**：当数据来自不同的子群体或存在不同的变异来源时，可能会出现过度离势。比如在调查不同学校学生的考试成绩时，学校之间的差异以及学生个体之间的差异可能导致成绩数据的离散程度较大，超过了基于单一分布假设下的预期。通过计算过度离势，可以帮助识别数据中的这种异质性，进而采取相应的分析方法，如分层分析或混合效应模型等。\n",
    "- **检验统计假设**：在一些统计检验中，如基于泊松分布或二项分布的假设检验，如果存在过度离势，可能会影响检验的准确性和有效性。计算过度离势可以帮助判断是否需要对检验方法进行调整，以获得更可靠的统计推断结果。例如，在对基因突变频率进行统计分析时，若存在过度离势，可能需要使用考虑了过度离势的修正检验方法，以避免得出错误的结论。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a414b962-dea1-497d-b729-52a83f68867a",
   "metadata": {},
   "source": [
    "在纵向数据分析中，重复测量数据是指对同一组个体在多个时间点上进行观测得到的数据。由于同一组个体的观测值之间存在内在的群聚特性，即个体自身的一些特质会使得其多次测量值之间具有相关性，这就可能导致数据出现过度离势的情况。以下是一个具体的例子：\n",
    "\n",
    "**研究背景**\n",
    "\n",
    "假设研究人员想要研究某种药物对高血压患者血压的影响。选取了50名高血压患者作为研究对象，在患者服用药物后的第1周、第2周、第3周、第4周分别测量他们的收缩压。\n",
    "\n",
    "**数据特征**\n",
    "\n",
    "- 一般来说，在没有其他因素干扰的情况下，按照常规的统计模型假设，血压值的波动应该在一个相对稳定的范围内，即数据的离散程度应该符合某种特定的分布（如正态分布）。\n",
    "- 然而，在这个研究中，每个患者自身的身体状况、生活习惯、对药物的吸收和反应等因素是相对稳定且具有个体特异性的，这些因素会导致同一个患者在不同时间点的血压测量值之间存在相关性。\n",
    "\n",
    "**导致过度离势的原因**\n",
    "\n",
    "- 例如，有的患者可能因为自身的遗传因素、生活方式（如长期坚持锻炼、饮食健康）等，使得其血压在整个观察期内都相对稳定且较低；而有的患者可能由于工作压力大、不按时服药等原因，血压波动较大且整体处于较高水平。\n",
    "- 这种个体间的差异会使得数据的离散程度大于基于常规统计模型（如假设数据独立同分布的模型）所预期的离散程度，即出现过度离势的现象。\n",
    "\n",
    "**数据表现**\n",
    "\n",
    "- 如果绘制出这些患者收缩压的箱线图，可以发现不同患者的血压数据分布差异较大，有些患者的血压值集中在较低的区间，而有些患者的血压值则集中在较高的区间，而且每个患者不同时间点的血压值也并非完全随机波动，而是呈现出一定的聚集性。\n",
    "- 从统计数据上看，计算得到的数据方差会明显大于假设数据独立时的理论方差，这就是数据存在过度离势的直观表现。\n",
    "\n",
    "在这个例子中，由于对同一组患者进行多次测量，患者个体的内在因素导致了数据的群聚特性，进而引发了过度离势。在分析这类数据时，如果不考虑过度离势，可能会导致模型拟合不准确，对药物效果的评估也会产生偏差。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "0e5d05fb-3258-407d-882f-bfa35f3dc927",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "10.1717029935069"
      ],
      "text/latex": [
       "10.1717029935069"
      ],
      "text/markdown": [
       "10.1717029935069"
      ],
      "text/plain": [
       "[1] 10.1717"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "deviance(fit)/df.residual(fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2aa3cdd1-9ae2-44d1-b55d-e0d11f826a2a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "package 'qcc' successfully unpacked and MD5 sums checked\n",
      "\n",
      "The downloaded binary packages are in\n",
      "\tC:\\Users\\xie.xiaokang\\AppData\\Local\\Temp\\RtmpQtKRr6\\downloaded_packages\n"
     ]
    }
   ],
   "source": [
    "# options(repos = c(CRAN = \"https://mirrors.aliyun.com/CRAN/\"))\n",
    "# install.packages(\"qcc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "f0ac8c3f-b5d4-4196-8536-def8dffacac2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "\"package 'qcc' was built under R version 4.4.3\"\n",
      "Package 'qcc' version 2.7\n",
      "\n",
      "Type 'citation(\"qcc\")' for citing this R package in publications.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                   \n",
      "Overdispersion test Obs.Var/Theor.Var Statistic p-value\n",
      "       poisson data          62.87013  3646.468       0\n"
     ]
    }
   ],
   "source": [
    "library(qcc)\n",
    "print(qcc.overdispersion.test(breslow.dat$sumY,type=\"poisson\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "6db7d695-3c9f-41af-ad40-e6a17c4b091e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "62.8701326259947"
      ],
      "text/latex": [
       "62.8701326259947"
      ],
      "text/markdown": [
       "62.8701326259947"
      ],
      "text/plain": [
       "[1] 62.87013"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lambda = mean(breslow.dat$sumY)\n",
    "variance = var(breslow.dat$sumY)\n",
    "variance/lambda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "fab4d423-82af-47a8-bbed-246ae81ed1ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = sumY ~ Base + Age + Trt, family = quasipoisson(), \n",
       "    data = breslow.dat)\n",
       "\n",
       "Coefficients:\n",
       "              Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept)   1.948826   0.465091   4.190 0.000102 ***\n",
       "Base          0.022652   0.001747  12.969  < 2e-16 ***\n",
       "Age           0.022740   0.013800   1.648 0.105085    \n",
       "Trtprogabide -0.152701   0.163943  -0.931 0.355702    \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for quasipoisson family taken to be 11.76075)\n",
       "\n",
       "    Null deviance: 2122.73  on 58  degrees of freedom\n",
       "Residual deviance:  559.44  on 55  degrees of freedom\n",
       "AIC: NA\n",
       "\n",
       "Number of Fisher Scoring iterations: 5\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fit.od <- glm(sumY~Base+Age+Trt,data=breslow.dat,family=quasipoisson())\n",
    "summary(fit.od)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a4f85d4-244a-4bd8-a77f-f9b4d3afc7bf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
